#Setup the working environment
setwd("/hpc/group/pbenfeylab/mz187/")
# Load the packages
suppressMessages(library(Seurat))
suppressMessages(library(SeuratWrappers))
suppressMessages(library(monocle3))
# Read rice root atlast Copilot project
rc.integrated <- readRDS("Rice_12S_WT_atlas_20221117_seu3.rds")
# Show the UMAP and cell number for each cell type cluster
DimPlot(rc.integrated, reduction = "umap",group.by = "celltype.anno")
table(rc.integrated$celltype.anno)
Stem Cell Niche Sclerenchyma Root Cap Atrichoblast Trichoblast
7746 8493 4724 4675 5401
Cortex Endodermis Phloem Pericycle Exodermis
17257 4362 2455 5088 10326
Protoxylem Metaxylem
7928 1288
# Define meristem, trichoblast and atrichoblast together for epidermis developmental trajectory analysis
rc.integrated$epidermis <- rc.integrated$celltype.anno %in% c("Stem Cell Niche","Atrichoblast","Trichoblast")
table(rc.integrated$epidermis)
table(rc.integrated$celltype.anno)
rc.integrated$epidermis <- as.character(rc.integrated$epidermis)
rc.integrated$epidermis <- gsub("TRUE", "Epidermal", rc.integrated$epidermis,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$epidermis <- gsub("FALSE", "Others", rc.integrated$epidermis,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$epidermis <- factor(rc.integrated$epidermis ,
levels=c("Epidermal", "Others"))
table(rc.integrated$epidermis, rc.integrated$celltype.anno)
# Get cell number included in each group
table(rc.integrated$epidermis)
FALSE TRUE 61921 17822
Stem Cell Niche Sclerenchyma Root Cap Atrichoblast Trichoblast
7746 8493 4724 4675 5401
Cortex Endodermis Phloem Pericycle Exodermis
17257 4362 2455 5088 10326
Protoxylem Metaxylem
7928 1288
Stem Cell Niche Sclerenchyma Root Cap Atrichoblast Trichoblast
Epidermal 7746 0 0 4675 5401
Others 0 8493 4724 0 0
Cortex Endodermis Phloem Pericycle Exodermis Protoxylem Metaxylem
Epidermal 0 0 0 0 0 0 0
Others 17257 4362 2455 5088 10326 7928 1288
Epidermal Others
17822 61921
## Subgroup the normal group for correlation analysis
rc.epidermis <- subset(rc.integrated, cells=colnames(rc.integrated)[which(rc.integrated$epidermis=="Epidermal")])
# Check the epidermis UMAP
options(repr.plot.width=8, repr.plot.height=6)
DimPlot(rc.epidermis, reduction = "umap",group.by = "celltype.anno")
# Rerun UMAP to make the clusters simpler
rc.epidermis <- RunUMAP(rc.epidermis, reduction = "pca", dims = 1:3)
rc.epidermis <- FindNeighbors(rc.epidermis, reduction = "pca", dims = 1:3)
rc.epidermis <- FindClusters(rc.epidermis, resolution = 0.5)
Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 16:01:05 UMAP embedding parameters a = 0.9922 b = 1.112 16:01:05 Read 17822 rows and found 3 numeric columns 16:01:06 Using Annoy for neighbor search, n_neighbors = 30 16:01:06 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 16:01:08 Writing NN index file to temp file /tmp/RtmpjuEspD/file20a4e066125b16 16:01:08 Searching Annoy index using 1 thread, search_k = 3000 16:01:16 Annoy recall = 100% 16:01:17 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 16:01:20 Initializing from normalized Laplacian + noise (using irlba) 16:01:23 Commencing optimization for 200 epochs, with 501752 positive edges 16:01:31 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 17822 Number of edges: 461212 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9335 Number of communities: 18 Elapsed time: 1 seconds
# Use Dimplot to check the clustering result
DimPlot(rc.epidermis, reduction = "umap",label = TRUE)
table(rc.epidermis$seurat_clusters)
DimPlot(rc.epidermis, reduction = "umap",group.by = "celltype.anno")
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 2004 1544 1427 1136 1112 1081 1046 1033 1032 1013 946 919 894 598 575 567 16 17 530 365
# Remove the "stem cell niche" cells which seems not relevant to epidermis
rc.epidermis$use <- rc.epidermis$seurat_clusters %in% c("1","3","4","8","12","15","17")
table(rc.epidermis$use)
table(rc.epidermis$seurat_clusters)
rc.epidermis$use <- as.character(rc.epidermis$use)
rc.epidermis$use <- gsub("TRUE", "nonuse", rc.epidermis$use,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.epidermis$use <- gsub("FALSE", "inuse", rc.epidermis$use,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.epidermis$use <- factor(rc.epidermis$use,
levels=c("inuse", "nonuse"))
table(rc.epidermis$use, rc.epidermis$seurat_clusters)
FALSE TRUE 11172 6650
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 2004 1544 1427 1136 1112 1081 1046 1033 1032 1013 946 919 894 598 575 567 16 17 530 365
0 1 2 3 4 5 6 7 8 9 10 11 12 13
inuse 2004 0 1427 0 0 1081 1046 1033 0 1013 946 919 0 598
nonuse 0 1544 0 1136 1112 0 0 0 1032 0 0 0 894 0
14 15 16 17
inuse 575 0 530 0
nonuse 0 567 0 365
## Subgroup the wanted stem cell nich, atrichoblast, trichoblast cells for monocle analysis
sub.epidermis <- subset(rc.epidermis, cells=colnames(rc.epidermis)[which(rc.epidermis$use=="inuse")])
# Rerun UMAP to make the clusters simpler
sub.epidermis <- RunUMAP(sub.epidermis, reduction = "pca", dims = 1:7)
sub.epidermis <- FindNeighbors(sub.epidermis, reduction = "pca", dims = 1:7)
sub.epidermis <- FindClusters(sub.epidermis, resolution = 0.1)
# Use Dimplot to check the clustering result
DimPlot(sub.epidermis, reduction = "umap",label = TRUE)
DimPlot(sub.epidermis, reduction = "umap",group.by = "celltype.anno")
16:02:53 UMAP embedding parameters a = 0.9922 b = 1.112 16:02:53 Read 11172 rows and found 7 numeric columns 16:02:53 Using Annoy for neighbor search, n_neighbors = 30 16:02:53 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 16:02:54 Writing NN index file to temp file /tmp/RtmpjuEspD/file20a4e03a9f2840 16:02:54 Searching Annoy index using 1 thread, search_k = 3000 16:02:58 Annoy recall = 100% 16:02:59 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 16:03:02 Initializing from normalized Laplacian + noise (using irlba) 16:03:02 Commencing optimization for 200 epochs, with 425024 positive edges 16:03:08 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 11172 Number of edges: 318177 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9623 Number of communities: 6 Elapsed time: 1 seconds
# There is still a separated Stem Cell Niche cluster, far from the epidermal cell lineages. Will do another subgroup to exclude them.
# Remove more "stem cell niche" cells which seems not relevant to epidermis
sub.epidermis$filter <- sub.epidermis$seurat_clusters %in% c("5")
table(sub.epidermis$filter)
table(sub.epidermis$seurat_clusters)
sub.epidermis$filter <- as.character(sub.epidermis$filter)
sub.epidermis$filter <- gsub("TRUE", "exclude", sub.epidermis$filter,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
sub.epidermis$filter <- gsub("FALSE", "include", sub.epidermis$filter,
ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
sub.epidermis$filter <- factor(sub.epidermis$filter,
levels=c("include", "exclude"))
table(sub.epidermis$filter, sub.epidermis$seurat_clusters)
FALSE TRUE 10642 530
0 1 2 3 4 5 2927 2015 1995 1867 1838 530
0 1 2 3 4 5
include 2927 2015 1995 1867 1838 0
exclude 0 0 0 0 0 530
## Further subgroup the wanted stem cell nich, atrichoblast, trichoblast cells for monocle analysis
sub.epidermis2 <- subset(sub.epidermis, cells=colnames(sub.epidermis)[which(sub.epidermis$filter=="include")])
# Read rice root subgrouped Copilot project
sub.epidermis2 <- readRDS("./1-Root_Atlas/sub.epidermis.rds")
# Re-run PCA for the data subet
sub.epidermis2 <- RunPCA(sub.epidermis2, npcs = 50, verbose = FALSE, approx = TRUE)
# Rerun UMAP to make the clusters simpler
sub.epidermis2 <- RunUMAP(sub.epidermis2, reduction = "pca",n.neighbors = 50L, min.dist = 0.4, dims = 1:7)
sub.epidermis2 <- FindNeighbors(sub.epidermis2, reduction = "pca", dims = 1:7)
sub.epidermis2 <- FindClusters(sub.epidermis2, resolution = 0.5)
# Use Dimplot to check the clustering result
options(repr.plot.width=18, repr.plot.height=6)
DimPlot(sub.epidermis2, reduction = "umap",label = TRUE,group.by = c("seurat_clusters","celltype.anno","time.anno"), ncol = 3)
14:15:24 UMAP embedding parameters a = 0.7669 b = 1.223 14:15:24 Read 10642 rows and found 7 numeric columns 14:15:24 Using Annoy for neighbor search, n_neighbors = 50 14:15:24 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 14:15:26 Writing NN index file to temp file /tmp/RtmpHIVa7I/file2f10485810e279 14:15:26 Searching Annoy index using 1 thread, search_k = 5000 14:15:33 Annoy recall = 100% 14:15:34 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 50 14:15:36 Initializing from normalized Laplacian + noise (using irlba) 14:15:37 Commencing optimization for 200 epochs, with 647838 positive edges 14:15:44 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 10642 Number of edges: 315339 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9164 Number of communities: 15 Elapsed time: 2 seconds
# Save the rds object the PCA calculated
saveRDS(sub.epidermis2, file = "./1-Root_Atlas/sub.epidermis-2.rds")
# Read rice root subgrouped Copilot project
sub.epidermis2 <- readRDS("./1-Root_Atlas/sub.epidermis-2.rds")
# Use Dimplot to check the clustering result
options(repr.plot.width=18, repr.plot.height=6)
DimPlot(sub.epidermis2, reduction = "umap",label = TRUE,group.by = c("seurat_clusters","celltype.anno","time.anno"), ncol = 3)
## Manually assign the clusters to different cell groups with cell type and developmental stage ananotation
sub.epidermis2$defined_clusters <- rep("Maturation Trichoblast", ncol(sub.epidermis2))
sub.epidermis2$defined_clusters[which((sub.epidermis2$seurat_clusters==10))] = "Meristem Epidermis"
sub.epidermis2$defined_clusters[which(sub.epidermis2$seurat_clusters==3)] = "Elongation Epidermis"
sub.epidermis2$defined_clusters[which((sub.epidermis2$seurat_clusters==11)|(sub.epidermis2$seurat_clusters==4)|(sub.epidermis2$seurat_clusters==5))] = "Elongation Atrichoblast"
sub.epidermis2$defined_clusters[which((sub.epidermis2$seurat_clusters==14)|(sub.epidermis2$seurat_clusters==8))] = "Elongation Trichoblast"
sub.epidermis2$defined_clusters[which((sub.epidermis2$seurat_clusters==9)|(sub.epidermis2$seurat_clusters==12)|(sub.epidermis2$seurat_clusters==13))] = "Maturation Atrichoblast"
# Check whether the assigning works
table(sub.epidermis2$defined_clusters)
options(repr.plot.width=8, repr.plot.height=6)
DimPlot(sub.epidermis2, reduction = "umap",group.by = "defined_clusters")
Elongation Atrichoblast Elongation Epidermis Elongation Trichoblast
2084 890 1020
Maturation Atrichoblast Maturation Trichoblast Meristem Epidermis
1547 4562 539
# Define the "toMonocle3" function, code by Che-Wei Hsu
ToMonocle3 <- function(seurat_object,
scale_all = FALSE,
assay = "SCT",
reduction_for_projection = "PCA",
UMAP_cluster_slot = NULL){
if(scale_all){
message("Getting residuals for all Seurat genes in chosen assay slot and placing in scale.data")
seurat_genes <- rownames(seurat_object[[assay]])
remaining_genes <- setdiff(seurat_genes, rownames(seurat_object[[assay]]@scale.data))
if(assay == "SCT"){
seurat_object <- Seurat::GetResidual(seurat_object, features = remaining_genes, assay = assay, umi.assay = "RNA")
} else {
seurat_object <- Seurat::ScaleData(seurat_object, features = rownames(seurat_object[[assay]]))
}
}
#We prep the seurat object by creating gene loadings for ALL genes in the Seurat scale.data slot. This is done to allow downstream monocle3 functions on gene_modules to work appropriately.
message("Projecting gene loadings for all Seurat genes in scale.data slot")
seurat_object <- Seurat::ProjectDim(seurat_object, reduction = reduction_for_projection, assay = assay)
##################
message("Initializing CDS object")
#Extract Seurat's log-transformed values
expression_matrix <- Seurat::GetAssayData(seurat_object, assay = assay, slot = "data")
#Extract Seurat meta_data
meta_data <- seurat_object@meta.data
#Extract gene names from Seurat object SCT slot to make CDS
seurat_genes <- data.frame(gene_short_name = rownames(seurat_object[[assay]]),
row.names = rownames(seurat_object[[assay]]))
new_cds <- monocle3::new_cell_data_set(expression_data = expression_matrix, cell_metadata = meta_data, gene_metadata = seurat_genes)
##################
message("Making an SCE object from the Seurat object to facilitate transfer of information from SCE to CDS")
sce <- as.SingleCellExperiment(seurat_object, assay = assay)
message("Loading in all Seurat reductions (PCA, HARMONY, UMAP, etc.) into CDS")
SingleCellExperiment::reducedDims(new_cds) <- SingleCellExperiment::reducedDims(sce)
message("Loading in specified Seurat assay into CDS")
SummarizedExperiment::assays(new_cds) <- SummarizedExperiment::assays(sce)
message("Loading in Seurat gene names into CDS")
SummarizedExperiment::rowData(new_cds) <- SummarizedExperiment::rowData(sce)
SummarizedExperiment::rowData(new_cds)$gene_short_name <- row.names(new_cds)
message("Loading in Seurat gene loadings into CDS")
new_cds@preprocess_aux$gene_loadings <- seurat_object@reductions[[reduction_for_projection]]@feature.loadings.projected
##################
message("Get user specified selected clusters (or active idents) from Seurat and load into CDS")
if(is.null(UMAP_cluster_slot)){
list_cluster <- Idents(seurat_object)
} else {
Idents(seurat_object) <- UMAP_cluster_slot
list_cluster <- Idents(seurat_object)
}
new_cds@clusters[["UMAP"]]$clusters <- list_cluster
#The next two commands are run in order to allow "order_cells" to be run in monocle3
rownames(new_cds@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
colnames(SingleCellExperiment::reducedDims(new_cds)[["UMAP"]]) <- NULL
##################
message("Setting all cells as belonging to one partition (multiple partitions not supported yet)")
recreate_partition <- c(rep(1, length(new_cds@colData@rownames)))
names(recreate_partition) <- new_cds@colData@rownames
recreate_partition <- as.factor(recreate_partition)
new_cds@clusters[["UMAP"]]$partitions <- recreate_partition
##################
message("Done")
new_cds
}
# Use monocle3 embedded function to convert seurat project to monocle3 project
mon3.cds <- ToMonocle3(sub.epidermis2,scale_all = FALSE,
assay = "SCT",
reduction_for_projection = "umap",
UMAP_cluster_slot = NULL)
Projecting gene loadings for all Seurat genes in scale.data slot Warning message in print.DimReduc(x = redeuc, dims = dims.print, nfeatures = nfeatures.print, : “Only 2 dimensions have been computed.”
UMAP_ 1 Positive: LOC-Os07g29390, LOC-Os05g39580, LOC-Os05g19600, LOC-Os07g44550, LOC-Os06g08490, LOC-Os03g49270, LOC-Os05g39560, LOC-Os05g02510, LOC-Os05g51119, LOC-Os01g64840 LOC-Os04g32550, LOC-Os05g50710, LOC-Os10g30330, LOC-Os01g54930, LOC-Os05g44930, LOC-Os10g35436, LOC-Os10g19180, LOC-Os03g61470, LOC-Os10g30340, LOC-Os06g20150 Negative: LOC-Os07g12830, LOC-Os06g08880, LOC-Os10g28520, LOC-Os03g55070, LOC-Os03g26450, LOC-Os07g16364, LOC-Os01g49120, LOC-Os04g49530, LOC-Os01g49219, LOC-Os03g46740 LOC-Os04g58700, LOC-Os01g49190, LOC-Os10g10720, LOC-Os10g32880, LOC-Os01g52810, LOC-Os03g43910, LOC-Os08g43930, LOC-Os06g23800, LOC-Os07g47290, LOC-Os01g49154 UMAP_ 2 Positive: LOC-Os03g26450, LOC-Os12g30540, LOC-Os01g49150, LOC-Os01g07310, LOC-Os04g53190, LOC-Os08g43540, LOC-Os10g27180, LOC-Os03g07340, LOC-Os05g15240, LOC-Os10g06580 LOC-Os01g49120, LOC-Os05g14740, LOC-Os02g44599, LOC-Os01g49154, LOC-Os04g25454, LOC-Os02g23939, LOC-Os09g30418, LOC-Os10g32880, LOC-Os06g08880, LOC-Os01g49210 Negative: LOC-Os11g24170, LOC-Os12g41710, LOC-Os12g41690, LOC-Os02g04730, LOC-Os10g05180, LOC-Os01g31220, LOC-Os03g36490, LOC-Os05g39560, LOC-Os02g42360, LOC-Os12g14059 LOC-Os07g28340, LOC-Os12g44020, LOC-Os02g06380, LOC-Os02g47970, LOC-Os02g41590, LOC-Os07g02210, LOC-Os05g02510, LOC-Os01g45020, LOC-Os12g35560, LOC-Os03g03610
Warning message: “Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.umap; see ?make.names for more details on syntax validity” Initializing CDS object Making an SCE object from the Seurat object to facilitate transfer of information from SCE to CDS Loading in all Seurat reductions (PCA, HARMONY, UMAP, etc.) into CDS Loading in specified Seurat assay into CDS Loading in Seurat gene names into CDS Loading in Seurat gene loadings into CDS Get user specified selected clusters (or active idents) from Seurat and load into CDS Setting all cells as belonging to one partition (multiple partitions not supported yet) Done
# Skip the cluster_cells process and run learn_graph function
mon3.cds <- learn_graph(mon3.cds)
# a helper function by Che-Wei to identify the root principal points:
get_earliest_principal_node <- function(cds, root="Meristem Epidermis"){
cell_ids <- which(colData(cds)[, "defined_clusters"] == root)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
mon3.cds <- order_cells(mon3.cds, root_pr_nodes=get_earliest_principal_node(mon3.cds))
|======================================================================| 100%
Warning message in igraph::graph.dfs(stree_ori, root = root_cell, neimode = "all", : “Argument `neimode' is deprecated; use `mode' instead”
# check the trajetory plot
options(repr.plot.width=8, repr.plot.height=8)
plot_cells(mon3.cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=TRUE,
graph_label_size=1.5,
cell_size = 1.5)
Warning message: “`select_()` was deprecated in dplyr 0.7.0. ℹ Please use `select()` instead. ℹ The deprecated feature was likely used in the monocle3 package. Please report the issue to the authors.”
## Conclusion: It seems that the Monocle3 can not identify the differentiation trajectory for trichoblast and the trichoblast and atrichoblast differentiation
## are connected, which is problematic. We thus decided to do the peusodotime analysis separately for trichoblast and atrichoblast, and then merge them.
## Monocle3 analysis, copied: https://stuartlab.org/signac/articles/monocle.html
epidermis.cds <- as.cell_data_set(sub.epidermis2)
Warning message: “Monocle 3 trajectories require cluster partitions, which Seurat does not calculate. Please run 'cluster_cells' on your cell_data_set object”
# More Monocle3 processing of epidermis cell dataset-2
epidermis.cds <- cluster_cells(cds = epidermis.cds, reduction_method = "UMAP", cluster_method = "louvain")
# More Monocle3 processing of epidermis cell dataset-3
epidermis.cds <- learn_graph(epidermis.cds, use_partition = TRUE)
|======================================================================| 100%
Warning message in igraph::graph.dfs(stree_ori, root = root_cell, neimode = "all", : “Argument `neimode' is deprecated; use `mode' instead”
Error: names(cds@clusters[[reduction_method]]$partitions) == NULL Traceback: 1. learn_graph(epidermis.cds, use_partition = TRUE) 2. project2MST(cds, project_point_to_line_segment, orthogonal_proj_tip, . verbose, reduction_method, rge_res_Y) 3. assertthat::assert_that(!is.null(names(cds@clusters[[reduction_method]]$partitions)), . msg = "names(cds@clusters[[reduction_method]]$partitions) == NULL")
# Monocle3-order cells
epidermis.cds <- order_cells(epidermis.cds, reduction_method = "UMAP", root_cells = colnames(epidermis.cds)[which(epidermis.cds$seurat_clusters == "13")])
## check the trajetory plot
options(repr.plot.width=8, repr.plot.height=6)
plot_cells(epidermis.cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=TRUE,
graph_label_size=1.5,
cell_size = 1.5)
# Extract the pseudotime values and add to the Seurat object
sub.epidermis2 <- AddMetaData(object = sub.epidermis2,
metadata = epidermis.cds@principal_graph_aux@listData$UMAP$pseudotime,
col.name = "pseudotimevalue")
# Plotting with converting UMAP 1 and UMAP2
options(repr.plot.width=6, repr.plot.height=6)
FeaturePlot(sub.epidermis2, "pseudotimevalue", pt.size = 0.3,dims = c(2, 1),cols = c("goldenrod1", "blue4"))
# Subgroup the atrichobladst cells
atri <- subset(sub.epidermis2, cells=colnames(sub.epidermis2)[which(sub.epidermis2$defined_clusters=="Meristem Epidermis"|
sub.epidermis2$defined_clusters=="Elongation Epidermis"|
sub.epidermis2$defined_clusters=="Elongation Atrichoblast"|
sub.epidermis2$defined_clusters=="Maturation Atrichoblast")])
# Check whether the subgrouping is correct
DimPlot(atri, reduction = "umap",group.by = "defined_clusters")
# Re-run PCA for the data subet
atri <- RunPCA(atri, npcs = 10, verbose = FALSE, approx = TRUE)
# Rerun UMAP to make the clusters simpler
atri <- RunUMAP(atri, reduction = "pca",n.neighbors = 100L, min.dist = 0.3, dims = 1:4)
atri <- FindNeighbors(atri, reduction = "pca", dims = 1:4)
atri <- FindClusters(atri, resolution = 0.5)
# Use Dimplot to check the clustering result
options(repr.plot.width=18, repr.plot.height=6)
DimPlot(atri, reduction = "umap",label = TRUE,group.by = c("seurat_clusters","defined_clusters","time.anno"), ncol = 3)
Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 15:31:26 UMAP embedding parameters a = 0.9922 b = 1.112 15:31:26 Read 5060 rows and found 4 numeric columns 15:31:26 Using Annoy for neighbor search, n_neighbors = 100 15:31:26 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 15:31:27 Writing NN index file to temp file /tmp/RtmpAVaWtp/file2fb33a7cb16696 15:31:27 Searching Annoy index using 1 thread, search_k = 10000 15:31:35 Annoy recall = 100% 15:31:37 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 100 15:31:41 Initializing from normalized Laplacian + noise (using irlba) 15:31:42 Commencing optimization for 500 epochs, with 430966 positive edges 15:31:56 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 5060 Number of edges: 137029 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9054 Number of communities: 12 Elapsed time: 0 seconds
# Use monocle3 embedded function to convert seurat project to monocle3 project
atri.cds <- ToMonocle3(atri,scale_all = FALSE,
assay = "SCT",
reduction_for_projection = "umap",
UMAP_cluster_slot = NULL)
Projecting gene loadings for all Seurat genes in scale.data slot Warning message in print.DimReduc(x = redeuc, dims = dims.print, nfeatures = nfeatures.print, : “Only 2 dimensions have been computed.”
UMAP_ 1 Positive: LOC-Os08g06344, LOC-Os04g38270, LOC-Os05g46430, LOC-Os08g01870, LOC-Os03g15070, LOC-Os01g05694, LOC-Os08g01600, LOC-Os02g33610, LOC-Os05g18790, LOC-Os04g25410 LOC-Os08g24310, LOC-Os01g67950, LOC-Os03g54780, LOC-Os05g51119, LOC-Os12g08820, LOC-Os08g29669, LOC-Os02g54120, LOC-Os03g02900, LOC-Os02g10900, LOC-Os09g36976 Negative: LOC-Os11g24170, LOC-Os03g51250, LOC-Os06g08620, LOC-Os05g02510, LOC-Os12g41690, LOC-Os03g43720, LOC-Os05g39560, LOC-Os03g03610, LOC-Os12g41710, LOC-Os01g55190 LOC-Os03g13140, LOC-Os02g04730, LOC-Os01g31220, LOC-Os07g28340, LOC-Os03g61990, LOC-Os12g14059, LOC-Os06g08490, LOC-Os03g15330, LOC-Os05g02520, LOC-Os02g43860 UMAP_ 2 Positive: LOC-Os07g29390, LOC-Os06g08490, LOC-Os07g44550, LOC-Os05g39580, LOC-Os05g19600, LOC-Os01g54930, LOC-Os01g64840, LOC-Os03g49270, LOC-Os10g19180, LOC-Os01g40470 LOC-Os10g30340, LOC-Os10g37600, LOC-Os10g30330, LOC-Os07g15640, LOC-Os05g27320, LOC-Os05g44930, LOC-Os01g45070, LOC-Os03g61470, LOC-Os10g35436, LOC-Os05g50710 Negative: LOC-Os12g41710, LOC-Os11g24170, LOC-Os12g41690, LOC-Os10g05180, LOC-Os12g14059, LOC-Os04g49530, LOC-Os05g39110, LOC-Os07g12830, LOC-Os01g31220, LOC-Os07g16364 LOC-Os02g47970, LOC-Os10g28520, LOC-Os01g55150, LOC-Os04g49610, LOC-Os01g02090, LOC-Os03g13160, LOC-Os02g42360, LOC-Os03g61990, LOC-Os07g02210, LOC-Os03g18110
Warning message: “Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.umap; see ?make.names for more details on syntax validity” Initializing CDS object Making an SCE object from the Seurat object to facilitate transfer of information from SCE to CDS Loading in all Seurat reductions (PCA, HARMONY, UMAP, etc.) into CDS Loading in specified Seurat assay into CDS Loading in Seurat gene names into CDS Loading in Seurat gene loadings into CDS Get user specified selected clusters (or active idents) from Seurat and load into CDS Setting all cells as belonging to one partition (multiple partitions not supported yet) Done
# Skip the cluster_cells process and run learn_graph function
atri.cds <- learn_graph(atri.cds)
# a helper function by Che-Wei to identify the root principal points:
get_earliest_principal_node <- function(cds, root="Meristem Epidermis"){
cell_ids <- which(colData(cds)[, "defined_clusters"] == root)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
atri.cds <- order_cells(atri.cds, root_pr_nodes=get_earliest_principal_node(atri.cds))
|======================================================================| 100%
Warning message in igraph::graph.dfs(stree_ori, root = root_cell, neimode = "all", : “Argument `neimode' is deprecated; use `mode' instead”
## check the trajetory plot
options(repr.plot.width=8, repr.plot.height=6)
plot_cells(atri.cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=TRUE,
graph_label_size=1.5,
cell_size = 1.5)
Warning message: “`select_()` was deprecated in dplyr 0.7.0. ℹ Please use `select()` instead. ℹ The deprecated feature was likely used in the monocle3 package. Please report the issue to the authors.”
# Subgroup the trichobladst cells
tri <- subset(sub.epidermis2, cells=colnames(sub.epidermis2)[which(sub.epidermis2$defined_clusters=="Meristem Epidermis"|
sub.epidermis2$defined_clusters=="Elongation Epidermis"|
sub.epidermis2$defined_clusters=="Elongation Trichoblast"|
sub.epidermis2$defined_clusters=="Maturation Trichoblast")])
# Check whether the subgrouping is correct
DimPlot(tri, reduction = "umap",group.by = "defined_clusters")
# Re-run PCA for the data subet
tri <- RunPCA(tri, npcs = 10, verbose = FALSE, approx = TRUE)
# Rerun UMAP to make the clusters simpler
tri <- RunUMAP(tri, reduction = "pca",n.neighbors = 100L, min.dist = 0.3, dims = 1:4)
tri <- FindNeighbors(tri, reduction = "pca", dims = 1:4)
tri <- FindClusters(tri, resolution = 0.5)
# Use Dimplot to check the clustering result
# Use Dimplot to check the clustering result
options(repr.plot.width=18, repr.plot.height=6)
DimPlot(tri, reduction = "umap",label = TRUE,group.by = c("seurat_clusters","defined_clusters","time.anno"), ncol = 3)
15:33:46 UMAP embedding parameters a = 0.9922 b = 1.112 15:33:46 Read 7011 rows and found 4 numeric columns 15:33:46 Using Annoy for neighbor search, n_neighbors = 100 15:33:46 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 15:33:48 Writing NN index file to temp file /tmp/RtmpAVaWtp/file2fb33a6dafeabf 15:33:48 Searching Annoy index using 1 thread, search_k = 10000 15:34:00 Annoy recall = 100% 15:34:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 100 15:34:06 Initializing from normalized Laplacian + noise (using irlba) 15:34:07 Commencing optimization for 500 epochs, with 576382 positive edges 15:34:25 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 7011 Number of edges: 190769 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9096 Number of communities: 15 Elapsed time: 1 seconds
# Use monocle3 embedded function to convert seurat project to monocle3 project
tri.cds <- ToMonocle3(tri,scale_all = FALSE,
assay = "SCT",
reduction_for_projection = "umap",
UMAP_cluster_slot = NULL)
Projecting gene loadings for all Seurat genes in scale.data slot Warning message in print.DimReduc(x = redeuc, dims = dims.print, nfeatures = nfeatures.print, : “Only 2 dimensions have been computed.”
UMAP_ 1 Positive: LOC-Os06g08880, LOC-Os03g26450, LOC-Os03g55070, LOC-Os01g49120, LOC-Os10g28520, LOC-Os07g12830, LOC-Os07g16364, LOC-Os04g49530, LOC-Os01g49219, LOC-Os03g46740 LOC-Os06g08620, LOC-Os01g49190, LOC-Os04g58700, LOC-Os01g52810, LOC-Os01g49154, LOC-Os03g51250, LOC-Os10g32880, LOC-Os10g10720, LOC-Os07g47290, LOC-Os04g48130 Negative: LOC-Os08g06344, LOC-Os04g38270, LOC-Os03g15070, LOC-Os01g05694, LOC-Os05g39580, LOC-Os05g19600, LOC-Os08g01870, LOC-Os05g51119, LOC-Os05g46430, LOC-Os02g33610 LOC-Os08g01600, LOC-Os05g18790, LOC-Os07g29390, LOC-Os01g67950, LOC-Os03g54780, LOC-Os10g38600, LOC-Os04g25410, LOC-Os04g32550, LOC-Os02g48280, LOC-Os02g35630 UMAP_ 2 Positive: LOC-Os12g41710, LOC-Os11g24170, LOC-Os02g42360, LOC-Os02g47970, LOC-Os05g39110, LOC-Os03g36490, LOC-Os12g41690, LOC-Os12g35560, LOC-Os07g28910, LOC-Os03g57545 LOC-Os04g49610, LOC-Os10g05180, LOC-Os12g44020, LOC-Os02g04730, LOC-Os01g31220, LOC-Os03g18110, LOC-Os02g41590, LOC-Os04g59510, LOC-Os06g25605, LOC-Os03g61990 Negative: LOC-Os05g02520, LOC-Os10g06580, LOC-Os05g15240, LOC-Os01g49150, LOC-Os03g48110, LOC-Os04g53190, LOC-Os02g44599, LOC-Os09g30410, LOC-Os01g40760, LOC-Os03g07340 LOC-Os12g30540, LOC-Os09g34847, LOC-Os05g14740, LOC-Os04g25454, LOC-Os05g45990, LOC-Os07g48880, LOC-Os01g55190, LOC-Os01g07310, LOC-Os05g42150, LOC-Os09g30418
Warning message: “Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.umap; see ?make.names for more details on syntax validity” Initializing CDS object Making an SCE object from the Seurat object to facilitate transfer of information from SCE to CDS Loading in all Seurat reductions (PCA, HARMONY, UMAP, etc.) into CDS Loading in specified Seurat assay into CDS Loading in Seurat gene names into CDS Loading in Seurat gene loadings into CDS Get user specified selected clusters (or active idents) from Seurat and load into CDS Setting all cells as belonging to one partition (multiple partitions not supported yet) Done
# Skip the cluster_cells process and run learn_graph function
tri.cds <- learn_graph(tri.cds)
# a helper function by Che-Wei to identify the root principal points:
get_earliest_principal_node <- function(cds, root="Meristem Epidermis"){
cell_ids <- which(colData(cds)[, "defined_clusters"] == root)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
tri.cds <- order_cells(tri.cds, root_pr_nodes=get_earliest_principal_node(tri.cds))
|======================================================================| 100%
Warning message in igraph::graph.dfs(stree_ori, root = root_cell, neimode = "all", : “Argument `neimode' is deprecated; use `mode' instead”
## check the trajetory plot
options(repr.plot.width=8, repr.plot.height=6)
plot_cells(tri.cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=TRUE,
graph_label_size=1.5,
cell_size = 1.5)
# Store the pseudotime values from atri and tri sub-objects back to the sub.epidermis object, for the common clusters, calculate the means: (A + B)/2
sub.epidermis2$pseudotime <- rep(0, ncol(sub.epidermis2))
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Maturation Atrichoblast")] = pseudotime(atri.cds)[which(atri$defined_clusters=="Maturation Atrichoblast")]
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Elongation Atrichoblast")] = pseudotime(atri.cds)[which(atri$defined_clusters=="Elongation Atrichoblast")]
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Maturation Trichoblast")] = pseudotime(tri.cds)[which(tri$defined_clusters=="Maturation Trichoblast")]
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Elongation Trichoblast")] = pseudotime(tri.cds)[which(tri$defined_clusters=="Elongation Trichoblast")]
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Elongation Epidermis")] = (((pseudotime(tri.cds)[which(tri$defined_clusters=="Elongation Epidermis")]) + (pseudotime(atri.cds)[which(atri$defined_clusters=="Elongation Epidermis")]))/2)
sub.epidermis2$pseudotime[which(sub.epidermis2$defined_clusters=="Meristem Epidermis")] =(((pseudotime(tri.cds)[which(tri$defined_clusters=="Meristem Epidermis")]) + (pseudotime(atri.cds)[which(atri$defined_clusters=="Meristem Epidermis")]))/2)
# change the UMAP2 value to invert the UMAP plot
sub.epidermis2[["umap"]]@cell.embeddings[,2] <- sub.epidermis2[["umap"]]@cell.embeddings[,2]*-1
# Plotting
options(repr.plot.width=6, repr.plot.height=6)
FeaturePlot(sub.epidermis2, "pseudotime", pt.size = 0.3,cols = c("blue4","goldenrod1"))
# Save the rds object with Pseudotime
saveRDS(sub.epidermis2, file = "./1-Root_Atlas/sub.epidermis-with_Pseudotime.rds")
# Read subgrouped epidermis cell Copilot object
sub.epidermis2 <- readRDS("./1-Root_Atlas/sub.epidermis-with_Pseudotime.rds")
# Convert the COPILOT object to Monocle3 cds, need to run "define 'ToMonocle3' function first"
cds <- ToMonocle3(sub.epidermis2,scale_all = FALSE,
assay = "SCT",
reduction_for_projection = "umap",
UMAP_cluster_slot = NULL)
Projecting gene loadings for all Seurat genes in scale.data slot Warning message in print.DimReduc(x = redeuc, dims = dims.print, nfeatures = nfeatures.print, : “Only 2 dimensions have been computed.”
UMAP_ 1 Positive: LOC-Os07g29390, LOC-Os05g39580, LOC-Os05g19600, LOC-Os07g44550, LOC-Os06g08490, LOC-Os03g49270, LOC-Os05g39560, LOC-Os05g02510, LOC-Os05g51119, LOC-Os01g64840 LOC-Os04g32550, LOC-Os05g50710, LOC-Os10g30330, LOC-Os01g54930, LOC-Os05g44930, LOC-Os10g35436, LOC-Os10g19180, LOC-Os03g61470, LOC-Os10g30340, LOC-Os06g20150 Negative: LOC-Os07g12830, LOC-Os06g08880, LOC-Os10g28520, LOC-Os03g55070, LOC-Os03g26450, LOC-Os07g16364, LOC-Os01g49120, LOC-Os04g49530, LOC-Os01g49219, LOC-Os03g46740 LOC-Os04g58700, LOC-Os01g49190, LOC-Os10g10720, LOC-Os10g32880, LOC-Os01g52810, LOC-Os03g43910, LOC-Os08g43930, LOC-Os06g23800, LOC-Os07g47290, LOC-Os01g49154 UMAP_ 2 Positive: LOC-Os11g24170, LOC-Os12g41710, LOC-Os12g41690, LOC-Os02g04730, LOC-Os10g05180, LOC-Os01g31220, LOC-Os03g36490, LOC-Os05g39560, LOC-Os02g42360, LOC-Os12g14059 LOC-Os07g28340, LOC-Os12g44020, LOC-Os02g06380, LOC-Os02g47970, LOC-Os02g41590, LOC-Os07g02210, LOC-Os05g02510, LOC-Os01g45020, LOC-Os12g35560, LOC-Os03g03610 Negative: LOC-Os03g26450, LOC-Os12g30540, LOC-Os01g49150, LOC-Os01g07310, LOC-Os04g53190, LOC-Os08g43540, LOC-Os10g27180, LOC-Os03g07340, LOC-Os05g15240, LOC-Os10g06580 LOC-Os01g49120, LOC-Os05g14740, LOC-Os02g44599, LOC-Os01g49154, LOC-Os04g25454, LOC-Os02g23939, LOC-Os09g30418, LOC-Os10g32880, LOC-Os06g08880, LOC-Os01g49210
Warning message: “Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.umap; see ?make.names for more details on syntax validity” Initializing CDS object Making an SCE object from the Seurat object to facilitate transfer of information from SCE to CDS Loading in all Seurat reductions (PCA, HARMONY, UMAP, etc.) into CDS Loading in specified Seurat assay into CDS Loading in Seurat gene names into CDS Loading in Seurat gene loadings into CDS Get user specified selected clusters (or active idents) from Seurat and load into CDS Setting all cells as belonging to one partition (multiple partitions not supported yet) Done
# Run learn_graph and order_cell functions in Monocle3 to create the metadata of pseuotime,
# which can be then replaces with our calculated pseudotime value
cds <- learn_graph(cds)
|======================================================================| 100%
Warning message in igraph::graph.dfs(stree_ori, root = root_cell, neimode = "all", : “Argument `neimode' is deprecated; use `mode' instead”
# Run order_cell function
# a helper function by Che-Wei to identify the root principal points:
get_earliest_principal_node <- function(cds, root="Meristem Epidermis"){
cell_ids <- which(colData(cds)[, "defined_clusters"] == root)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
# Replace the Monocle3 pseudotime values with our calculated ones (calculate separately for trichoblast and atrichoblast)
cds@principal_graph_aux$UMAP$pseudotime <- sub.epidermis2$pseudotime
# based on the pseudotime value to define 10 groups
pseu <- pseudotime(cds)
s10 <- c(seq(1,round(ncol(cds)/10)*10,round(length(pseu)/10)),ncol(cds)+1)
pseudotime.group <- rep("NA", ncol(cds))
pseudotime.cell.name <- rep("NA", ncol(cds))
for (i in 1:10){
pseudotime.group[s10[i]: s10[i+1]-1] <- paste0("G",i-1)
pseudotime.cell.name[s10[i]: s10[i+1]-1] <- names(sort(rank(pseu))[s10[i]: s10[i+1]-1])
}
names(pseudotime.group) <- pseudotime.cell.name
part <- pseudotime.group
cell_group_df <- tibble::tibble(cell=row.names(colData(cds)),
cell_group=part[colnames(cds)])
cds$pseudotime.group <- cell_group_df$cell_group
suppressMessages(library(ggplot2))
suppressMessages(library(RColorBrewer))
brewer.pal(10,"Spectral")
# Plot the cell differentiation trajectory
options(repr.plot.width=8, repr.plot.height=8)
plot_cells(cds,
color_cells_by = "pseudotime.group",
show_trajectory_graph = FALSE,
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5,
cell_size = 1.5) + scale_colour_manual(values = brewer.pal(10,"Spectral"))
## Specify one epdiermis marker gene and three trichoblast marker gene
target_genes <- c("LOC-Os03g19990", "LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
# Check marker gene expression with featureplots
DefaultAssay(sub.epidermis2) <- "SCT"
options(repr.plot.width=12, repr.plot.height=12)
FeaturePlot(sub.epidermis2, features = target_genes,
cols = c("grey", "red"), label=F, repel=F, pt.size = 1,
label.size = 1.5, order=T, ncol = 2)
# Plot the gene expression dyanmics with monocle featureplot
options(repr.plot.width=12, repr.plot.height=12)
plot_cells(cds, genes=target_genes,
show_trajectory_graph=FALSE,
label_cell_groups=FALSE,
label_leaves=FALSE,cell_size = 1)
# Plot the gene expression dyanmics with linear plot
cds$pseudotime.group <- factor(cds$pseudotime.group, levels = rev(c("G0","G1","G2","G3","G4","G5","G6","G7","G8","G9")))
lineage_cds <- cds[rowData(cds)$gene_short_name %in% target_genes,]
options(repr.plot.width=12, repr.plot.height=6)
plot_genes_in_pseudotime(lineage_cds,
color_cells_by="pseudotime.group",
min_expr=0.5, ncol=2, panel_order = target_genes
)+ scale_colour_manual(values = rev(brewer.pal(10,"Spectral")))
# Identify the DE genes with Monocle3 graph_test function
test_res <- graph_test(
cds,
neighbor_graph = c("knn", "principal_graph"),
reduction_method = "UMAP",
k = 25,
method = c("Moran_I"),
alternative = "greater",
expression_family = "quasipoisson",
cores = 4,
verbose = FALSE
)
|====================================================| 100%, Elapsed 01:27:29
# Save the DE gene list
write.csv(test_res, "./1-Root_Atlas/test_res-by-Mond3.csv",row.names = TRUE)
# Import the DE gene list
test_res <- read.table(file = "./1-Root_Atlas/test_res-by-Mond3.csv",row.names = 1, header=TRUE, sep = ",")
# Heatmap plot to reveal the gene functions
deg_ids <- rownames(subset(test_res, q_value < 0.01))
deg_name <- subset(test_res, q_value < 0.01)$gene_short_name
head (test_res)
| status | p_value | morans_test_statistic | morans_I | gene_short_name | q_value | |
|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <chr> | <dbl> | |
| LOC-Os01g01010 | OK | 5.116516e-15 | 7.736327 | 0.020236653 | LOC-Os01g01010 | 1.230320e-14 |
| LOC-Os01g01019 | OK | 5.164650e-02 | 1.629095 | 0.004159613 | LOC-Os01g01019 | 7.416631e-02 |
| LOC-Os01g01030 | OK | 2.663745e-03 | 2.786535 | 0.007164346 | LOC-Os01g01030 | 4.262551e-03 |
| LOC-Os01g01040 | OK | 6.365610e-43 | 13.683555 | 0.035926369 | LOC-Os01g01040 | 2.203258e-42 |
| LOC-Os01g01050 | OK | 1.118794e-59 | 16.249892 | 0.042637319 | LOC-Os01g01050 | 4.375333e-59 |
| LOC-Os01g01060 | OK | 0.000000e+00 | 124.424655 | 0.327407121 | LOC-Os01g01060 | 0.000000e+00 |
head (deg_name)
## Pick the gene to be labeled, not working
genes_to_label <- filter(test_res, deg_name %in% c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050"))
# Use Monocle3 find_gene_modules function
gene_module_df <- find_gene_modules(
cds[deg_ids,],
reduction_method = c("UMAP"),
max_components = 2,
umap.metric = "cosine",
umap.min_dist = 0.1,
umap.n_neighbors = 15L,
umap.fast_sgd = FALSE,
umap.nn_method = "annoy",
k = 20,
leiden_iter = 1,
partition_qval = 0.05,
weight = FALSE,
resolution = 10^-6,
random_seed = 0L,
cores = 1,
verbose = FALSE
)
Spectral initialization failed to converge, using random initialization instead
# Heatmap plot to reveal the gene functions
draw_colnames_horizontal <- function (coln, gaps, ...) {
coord <- pheatmap:::find_coordinates(length(coln), gaps)
x <- coord$coord - 0.5 * coord$size
res <- grid::textGrob(
coln, x = x, y = unit(1, "npc") - unit(3,"bigpts"),
vjust = 1.5, hjust = 0.5, rot = 0, gp = grid::gpar(...)
)
return(res)
}
assignInNamespace(
x = "draw_colnames",
value = "draw_colnames_horizontal",
ns = asNamespace("pheatmap")
)
agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
colnames(agg_mat) <- stringr::str_c("", colnames(agg_mat))
options(repr.plot.width=10, repr.plot.height=10)
setHook("grid.newpage", function() pushViewport(viewport(x=1,y=1,width=0.9, height=0.9, name="vp", just=c("right","top"))), action="prepend")
pheatmap::pheatmap(agg_mat, cluster_rows=TRUE, cluster_cols=FALSE,
scale="column", clustering_method="complete",
fontsize=12, labels_row=rep("", length(unique(gene_module_df$module))))
setHook("grid.newpage", NULL, "replace")
grid.text("pseudotime group", y=-0.02, gp=grid::gpar(fontsize=16))
grid.text("gene module", x=-0.05, rot=90, gp=grid::gpar(fontsize=16))
options(repr.plot.width=10, repr.plot.height=50)
setHook("grid.newpage", function() pushViewport(viewport(x=1,y=1,width=0.9, height=0.9, name="vp", just=c("right","top"))), action="prepend")
pheatmap::pheatmap(agg_mat, cluster_rows=TRUE, cluster_cols=FALSE,
scale="column", clustering_method="complete",
fontsize=12)
setHook("grid.newpage", NULL, "replace")
grid.text("pseudotime group", y=-0.02, gp=grid::gpar(fontsize=16))
grid.text("gene module", x=-0.05, rot=90, gp=grid::gpar(fontsize=16))
library(ComplexHeatmap)
class(agg_mat)
show.label <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
show.id <- c("Module 40", "Module 11", "Module 99")
mark_anno = rowAnnotation(foo = anno_mark(at = match(show.id,rownames(agg_mat)), labels = show.label, labels_gp = gpar(col = "black", fontsize = 16)))
H <- Heatmap(as.matrix(agg_mat),
col = colorRampPalette(rev(brewer.pal(n = 11, name =
"RdYlBu")))(11),
show_column_names = F,
cluster_columns = F,
cluster_rows=T,
show_row_names = F,
right_annotation=mark_anno,
show_heatmap_legend=F,
clustering_method_rows="complete",
use_raster=F,
heatmap_legend_param = list(title_position="topcenter", title = "Expression", direction="horizontal"))
options(repr.plot.width=10, repr.plot.height=10)
draw(H)
show.label <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
show.id <- c("Module 40", "Module 11", "Module 99")
mark_anno = rowAnnotation(foo = anno_mark(at = match(show.id,rownames(agg_mat)), labels = show.label, labels_gp = gpar(col = "black", fontsize = 16)))
H <- Heatmap(scale(as.matrix(agg_mat)),
col = colorRampPalette(rev(brewer.pal(n = 11, name =
"RdYlBu")))(11),
show_column_names = F,
cluster_columns = F,
cluster_rows=T,
show_row_names = T,
right_annotation=mark_anno,
show_heatmap_legend=F,
clustering_method_rows="complete",
use_raster=F,
heatmap_legend_param = list(title_position="topcenter", title = "Expression", direction="horizontal"))
order <- rownames(agg_mat)[row_order(H)]
Warning message: “The heatmap has not been initialized. You might have different results if you repeatedly execute this function, e.g. when row_km/column_km was set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.”
match('Module 36', order)
new_order <- c(rev(order[1:58]),order[59:244])
agg_mat <- agg_mat[new_order,]
show.label <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
show.id <- c("Module 40", "Module 11", "Module 99")
mark_anno = rowAnnotation(foo = anno_mark(at = match(show.id,rownames(agg_mat)), labels = show.label, labels_gp = gpar(col = "black", fontsize = 16)))
H <- Heatmap(scale(as.matrix(agg_mat)),
col = colorRampPalette(rev(brewer.pal(n = 11, name =
"RdYlBu")))(11),
show_column_names = F,
cluster_columns = F,
cluster_rows=F,
show_row_names = F,
right_annotation=mark_anno,
show_heatmap_legend=F,
clustering_method_rows="complete",
use_raster=F,
heatmap_legend_param = list(title_position="topcenter", title = "Expression", direction="horizontal"))
options(repr.plot.width=10, repr.plot.height=10)
draw(H)
head(agg_mat)
[[ suppressing 10 column names ‘G0’, ‘G1’, ‘G2’ ... ]]
6 x 10 sparse Matrix of class "dgCMatrix"
Module 1 -0.1334537 -0.14754157 -0.3150538 -0.1767659 -0.07507645 -0.06467763
Module 2 0.2936735 -0.10560607 -0.3450147 -0.1036743 -0.10002747 -0.10695148
Module 3 -0.7208790 -0.69226044 -0.6460777 -0.2687351 -0.08150884 0.04492619
Module 4 0.2888668 -0.22825689 -0.4183544 -0.2137301 -0.14714879 -0.14648676
Module 5 0.2155738 0.26105565 -0.1678768 -0.1129208 -0.08103566 -0.20374199
Module 6 0.0003125 0.03490679 -0.2031616 -0.2008526 -0.08843608 -0.11795392
Module 1 0.09255397 0.27402839 0.25791969 0.27657421
Module 2 -0.04251751 0.20883607 0.15246992 0.14687426
Module 3 0.39267546 0.64025783 0.76471404 0.56461802
Module 4 0.05865392 0.27764798 0.32958883 0.19516682
Module 5 -0.05882884 0.01707834 0.03861445 0.08444121
Module 6 0.04469907 0.19795312 0.15513694 0.17292993
order <- names(sort(apply(scale(agg_mat), MARGIN = 1, which.max)))
agg_mat <- scale(agg_mat)[order,]
options(repr.plot.width=10, repr.plot.height=30)
setHook("grid.newpage", function() pushViewport(viewport(x=1,y=1,width=0.9, height=0.9, name="vp", just=c("right","top"))), action="prepend")
pheatmap::pheatmap(agg_mat, cluster_rows=FALSE, cluster_cols=FALSE,
scale="column", clustering_method="ward.D2",
fontsize=12)
setHook("grid.newpage", NULL, "replace")
grid.text("pseudotime group", y=-0.02, gp=grid::gpar(fontsize=16))
grid.text("gene module", x=-0.05, rot=90, gp=grid::gpar(fontsize=16))
## Specify one epdiermis marker gene and three trichoblast marker gene
target_genes <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
for (i in target_genes){
print(gene_module_df$module[which(gene_module_df$id==i)])
}
[1] 40 244 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 244 [1] 11 244 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 244 [1] 99 244 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 244
which(order == "Module 40")
which(order == "Module 11")
which(order == "Module 99")
head(gene_module_df)
| id | module | supermodule | dim_1 | dim_2 |
|---|---|---|---|---|
| <chr> | <fct> | <fct> | <dbl> | <dbl> |
| LOC-Os01g01010 | 182 | 1 | 5.350191 | 12.300430 |
| LOC-Os01g01030 | 67 | 1 | -10.561002 | 1.902842 |
| LOC-Os01g01040 | 75 | 1 | -6.356702 | -7.412914 |
| LOC-Os01g01050 | 16 | 1 | -3.345216 | -2.505375 |
| LOC-Os01g01060 | 203 | 1 | -2.348689 | 3.785970 |
| LOC-Os01g01070 | 10 | 1 | 7.493682 | 4.606720 |
head(cell_group_df)
| cell | cell_group |
|---|---|
| <chr> | <chr> |
| AAACCCAAGTTGCTCA_1 | G7 |
| AAACCCAGTACGACTT_1 | G7 |
| AAACGAAAGCGTGCTC_1 | G0 |
| AAACGAACACCCAACG_1 | G9 |
| AAACGAAGTGTCCATA_1 | G7 |
| AAACGAATCTTCTGGC_1 | G4 |
# Load ComplexHeatmap library
suppressMessages(library(ComplexHeatmap))
# Define a function to label specific genes on the heatmap
label_genes <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
draw_colnames_with_labels <- function (coln, gaps, label_genes, ...) {
coord <- pheatmap:::find_coordinates(length(coln), gaps)
x <- coord$coord - 0.5 * coord$size
res <- grid::textGrob(
coln, x = x, y = unit(1, "npc") - unit(3, "bigpts"),
vjust = 1.5, hjust = 0.5, rot = 0, gp = grid::gpar(...)
)
# Label specific genes on the right side
label_indices <- which(coln %in% label_genes)
res_label <- grid::textGrob(
label_genes, x = rep(max(x) + 1, length(label_genes)),
y = unit(1, "npc") - unit(3, "bigpts") - coord$size[label_indices] * seq_along(label_genes),
vjust = 1.5, hjust = 0, rot = 0, gp = grid::gpar(...)
)
return(grid::gList(res, res_label))
}
# Assign the new function to the pheatmap namespace
assignInNamespace(
x = "draw_colnames",
value = "draw_colnames_with_labels",
ns = asNamespace("pheatmap")
)
# Generate the heatmap with labeled genes
pheatmap::pheatmap(agg_mat, cluster_rows = TRUE, cluster_cols = FALSE,
scale = "column", clustering_method = "ward.D2",
fontsize = 12, labels_row = rep("", length(unique(gene_module_df$module))),
draw_colnames = TRUE, label_genes = label_genes)
# Additional formatting
setHook("grid.newpage", NULL, "replace")
grid.text("pseudotime group", y = -0.02, gp = grid::gpar(fontsize = 16))
grid.text("gene module", x = -0.05, rot = 90, gp = grid::gpar(fontsize = 16))
Error in as.unit(e2): object is not coercible to a unit
Traceback:
1. pheatmap::pheatmap(agg_mat, cluster_rows = TRUE, cluster_cols = FALSE,
. scale = "column", clustering_method = "ward.D2", fontsize = 12,
. labels_row = rep("", length(unique(gene_module_df$module))),
. draw_colnames = TRUE, label_genes = label_genes)
2. heatmap_motor(mat, border_color = border_color, cellwidth = cellwidth,
. cellheight = cellheight, treeheight_col = treeheight_col,
. treeheight_row = treeheight_row, tree_col = tree_col, tree_row = tree_row,
. filename = filename, width = width, height = height, breaks = breaks,
. color = color, legend = legend, annotation_row = annotation_row,
. annotation_col = annotation_col, annotation_colors = annotation_colors,
. annotation_legend = annotation_legend, annotation_names_row = annotation_names_row,
. annotation_names_col = annotation_names_col, main = main,
. fontsize = fontsize, fontsize_row = fontsize_row, fontsize_col = fontsize_col,
. hjust_col = hjust_col, vjust_col = vjust_col, angle_col = angle_col,
. fmat = fmat, fontsize_number = fontsize_number, number_color = number_color,
. gaps_row = gaps_row, gaps_col = gaps_col, labels_row = labels_row,
. labels_col = labels_col, ...)
3. do.call(draw_colnames, pars)
4. draw_colnames_with_labels(c("G0", "G1", "G2", "G3", "G4", "G5",
. "G6", "G7", "G8", "G9"), gaps = NULL, fontsize = 12, hjust_col = 0,
. vjust_col = 0.5, angle_col = 270, draw_colnames = TRUE, label_genes = c("LOC-Os12g05380",
. "LOC-Os10g42750", "LOC-Os06g48050"))
5. grid::textGrob(label_genes, x = rep(max(x) + 1, length(label_genes)),
. y = unit(1, "npc") - unit(3, "bigpts") - coord$size[label_indices] *
. seq_along(label_genes), vjust = 1.5, hjust = 0, rot = 0,
. gp = grid::gpar(...)) # at line 14-18 of file <text>
6. is.unit(x)
7. Ops.unit(max(x), 1) # at line 14-18 of file <text>
8. as.unit(e2)
# Genes to label on the right side
label_genes <- c("LOC-Os12g05380", "LOC-Os10g42750", "LOC-Os06g48050")
# Convert agg_mat to a data frame
agg_mat_df <- as.data.frame(agg_mat)
agg_mat_df$gene <- rownames(agg_mat_df)
# Create a ggplot object
heatmap_plot <- ggplot(data = agg_mat_df, aes(x = colnames(agg_mat_df), y = gene)) +
geom_tile(aes(fill = agg_mat_df[[1]]), color = "white") + # Assuming the first column represents values
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
axis.text.y = element_text(size = 8),
legend.position = "none") +
labs(title = "Gene Expression Heatmap",
x = "Gene",
y = "Module")
# Label specific genes on the right side
label_df <- data.frame(x = rep(length(colnames(agg_mat_df)) + 1, length(label_genes)),
y = seq_along(label_genes),
label = label_genes)
heatmap_plot +
geom_text(data = label_df, aes(x = x, y = y, label = label), vjust = 1, hjust = 0) +
theme(plot.margin = margin(0.5, 0.5, 0.5, 5, "cm")) # Adjust the right margin for gene labels
ERROR while rich displaying an object: Error in `geom_tile()`: ! Problem while computing aesthetics. ℹ Error occurred in the 1st layer. Caused by error in `check_aesthetics()`: ! Aesthetics must be either length 1 or the same as the data (245) ✖ Fix the following mappings: `x` Traceback: 1. tryCatch(withCallingHandlers({ . if (!mime %in% names(repr::mime2repr)) . stop("No repr_* for mimetype ", mime, " in repr::mime2repr") . rpr <- repr::mime2repr[[mime]](obj) . if (is.null(rpr)) . return(NULL) . prepare_content(is.raw(rpr), rpr) . }, error = error_handler), error = outer_handler) 2. tryCatchList(expr, classes, parentenv, handlers) 3. tryCatchOne(expr, names, parentenv, handlers[[1L]]) 4. doTryCatch(return(expr), name, parentenv, handler) 5. withCallingHandlers({ . if (!mime %in% names(repr::mime2repr)) . stop("No repr_* for mimetype ", mime, " in repr::mime2repr") . rpr <- repr::mime2repr[[mime]](obj) . if (is.null(rpr)) . return(NULL) . prepare_content(is.raw(rpr), rpr) . }, error = error_handler) 6. repr::mime2repr[[mime]](obj) 7. repr_text.default(obj) 8. paste(capture.output(print(obj)), collapse = "\n") 9. capture.output(print(obj)) 10. withVisible(...elt(i)) 11. print(obj) 12. print.ggplot(obj) 13. ggplot_build(x) 14. ggplot_build.ggplot(x) 15. by_layer(function(l, d) l$compute_aesthetics(d, plot), layers, . data, "computing aesthetics") 16. try_fetch(for (i in seq_along(data)) { . out[[i]] <- f(l = layers[[i]], d = data[[i]]) . }, error = function(cnd) { . cli::cli_abort(c("Problem while {step}.", i = "Error occurred in the {ordinal(i)} layer."), . call = layers[[i]]$constructor, parent = cnd) . }) 17. tryCatch(withCallingHandlers(expr, condition = function(cnd) { . { . .__handler_frame__. <- TRUE . .__setup_frame__. <- frame . if (inherits(cnd, "message")) { . except <- c("warning", "error") . } . else if (inherits(cnd, "warning")) { . except <- "error" . } . else { . except <- "" . } . } . while (!is_null(cnd)) { . if (inherits(cnd, "error")) { . out <- handlers[[1L]](cnd) . if (!inherits(out, "rlang_zap")) . throw(out) . } . inherit <- .subset2(.subset2(cnd, "rlang"), "inherit") . if (is_false(inherit)) { . return() . } . cnd <- .subset2(cnd, "parent") . } . }), stackOverflowError = handlers[[1L]]) 18. tryCatchList(expr, classes, parentenv, handlers) 19. tryCatchOne(expr, names, parentenv, handlers[[1L]]) 20. doTryCatch(return(expr), name, parentenv, handler) 21. withCallingHandlers(expr, condition = function(cnd) { . { . .__handler_frame__. <- TRUE . .__setup_frame__. <- frame . if (inherits(cnd, "message")) { . except <- c("warning", "error") . } . else if (inherits(cnd, "warning")) { . except <- "error" . } . else { . except <- "" . } . } . while (!is_null(cnd)) { . if (inherits(cnd, "error")) { . out <- handlers[[1L]](cnd) . if (!inherits(out, "rlang_zap")) . throw(out) . } . inherit <- .subset2(.subset2(cnd, "rlang"), "inherit") . if (is_false(inherit)) { . return() . } . cnd <- .subset2(cnd, "parent") . } . }) 22. f(l = layers[[i]], d = data[[i]]) 23. l$compute_aesthetics(d, plot) 24. compute_aesthetics(..., self = self) 25. check_aesthetics(evaled, n) 26. cli::cli_abort(c("Aesthetics must be either length 1 or the same as the data ({n})", . x = "Fix the following mappings: {.col {names(which(!good))}}")) 27. rlang::abort(message, ..., call = call, use_cli_format = TRUE, . .frame = .frame) 28. signal_abort(cnd, .file) 29. signalCondition(cnd) 30. (function (cnd) . { . { . .__handler_frame__. <- TRUE . .__setup_frame__. <- frame . if (inherits(cnd, "message")) { . except <- c("warning", "error") . } . else if (inherits(cnd, "warning")) { . except <- "error" . } . else { . except <- "" . } . } . while (!is_null(cnd)) { . if (inherits(cnd, "error")) { . out <- handlers[[1L]](cnd) . if (!inherits(out, "rlang_zap")) . throw(out) . } . inherit <- .subset2(.subset2(cnd, "rlang"), "inherit") . if (is_false(inherit)) { . return() . } . cnd <- .subset2(cnd, "parent") . } . })(structure(list(message = structure("Aesthetics must be either length 1 or the same as the data (245)", names = ""), . trace = structure(list(call = list(IRkernel::main(), kernel$run(), . handle_shell(), executor$execute(msg), tryCatch(evaluate(request$content$code, . envir = .GlobalEnv, output_handler = oh, stop_on_error = 1L), . interrupt = function(cond) { . log_debug("Interrupt during execution") . interrupted <<- TRUE . }, error = .self$handle_error), tryCatchList(expr, . classes, parentenv, handlers), tryCatchOne(tryCatchList(expr, . names[-nh], parentenv, handlers[-nh]), names[nh], . parentenv, handlers[[nh]]), doTryCatch(return(expr), . name, parentenv, handler), tryCatchList(expr, names[-nh], . parentenv, handlers[-nh]), tryCatchOne(expr, names, . parentenv, handlers[[1L]]), doTryCatch(return(expr), . name, parentenv, handler), evaluate(request$content$code, . envir = .GlobalEnv, output_handler = oh, stop_on_error = 1L), . evaluate_call(expr, parsed$src[[i]], envir = envir, enclos = enclos, . debug = debug, last = i == length(out), use_try = stop_on_error != . 2L, keep_warning = keep_warning, keep_message = keep_message, . log_echo = log_echo, log_warning = log_warning, output_handler = output_handler, . include_timing = include_timing), handle(pv <- withCallingHandlers(withVisible(value_fun(ev$value, . ev$visible)), warning = wHandler, error = eHandler, . message = mHandler)), try(f, silent = TRUE), tryCatch(expr, . error = function(e) { . call <- conditionCall(e) . if (!is.null(call)) { . if (identical(call[[1L]], quote(doTryCatch))) . call <- sys.call(-4L) . dcall <- deparse(call, nlines = 1L) . prefix <- paste("Error in", dcall, ": ") . LONG <- 75L . sm <- strsplit(conditionMessage(e), "\n")[[1L]] . w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], . type = "w") . if (is.na(w)) . w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], . type = "b") . if (w > LONG) . prefix <- paste0(prefix, "\n ") . } . else prefix <- "Error : " . msg <- paste0(prefix, conditionMessage(e), "\n") . .Internal(seterrmessage(msg[1L])) . if (!silent && isTRUE(getOption("show.error.messages"))) { . cat(msg, file = outFile) . .Internal(printDeferredWarnings()) . } . invisible(structure(msg, class = "try-error", . condition = e)) . }), tryCatchList(expr, classes, parentenv, handlers), . tryCatchOne(expr, names, parentenv, handlers[[1L]]), . doTryCatch(return(expr), name, parentenv, handler), withCallingHandlers(withVisible(value_fun(ev$value, . ev$visible)), warning = wHandler, error = eHandler, . message = mHandler), withVisible(value_fun(ev$value, . ev$visible)), value_fun(ev$value, ev$visible), prepare_mimebundle_kernel(obj, . .self$handle_display_error), prepare_mimebundle(obj, . "text/plain", error_handler = handle_display_error), . filter_map(mimetypes, function(mime) { . tryCatch(withCallingHandlers({ . if (!mime %in% names(repr::mime2repr)) . stop("No repr_* for mimetype ", mime, " in repr::mime2repr") . rpr <- repr::mime2repr[[mime]](obj) . if (is.null(rpr)) . return(NULL) . prepare_content(is.raw(rpr), rpr) . }, error = error_handler), error = outer_handler) . }), Filter(Negate(is.null), sapply(x, f, simplify = simplify)), . unlist(lapply(x, f)), lapply(x, f), sapply(x, f, simplify = simplify), . lapply(X = X, FUN = FUN, ...), FUN(X[[i]], ...), tryCatch(withCallingHandlers({ . if (!mime %in% names(repr::mime2repr)) . stop("No repr_* for mimetype ", mime, " in repr::mime2repr") . rpr <- repr::mime2repr[[mime]](obj) . if (is.null(rpr)) . return(NULL) . prepare_content(is.raw(rpr), rpr) . }, error = error_handler), error = outer_handler), tryCatchList(expr, . classes, parentenv, handlers), tryCatchOne(expr, . names, parentenv, handlers[[1L]]), doTryCatch(return(expr), . name, parentenv, handler), withCallingHandlers({ . if (!mime %in% names(repr::mime2repr)) . stop("No repr_* for mimetype ", mime, " in repr::mime2repr") . rpr <- repr::mime2repr[[mime]](obj) . if (is.null(rpr)) . return(NULL) . prepare_content(is.raw(rpr), rpr) . }, error = error_handler), repr::mime2repr[[mime]](obj), . repr_text.default(obj), paste(capture.output(print(obj)), . collapse = "\n"), capture.output(print(obj)), withVisible(...elt(i)), . print(obj), print.ggplot(obj), ggplot_build(x), ggplot_build.ggplot(x), . by_layer(function(l, d) l$compute_aesthetics(d, plot), . layers, data, "computing aesthetics"), try_fetch(for (i in seq_along(data)) { . out[[i]] <- f(l = layers[[i]], d = data[[i]]) . }, error = function(cnd) { . cli::cli_abort(c("Problem while {step}.", i = "Error occurred in the {ordinal(i)} layer."), . call = layers[[i]]$constructor, parent = cnd) . }), tryCatch(withCallingHandlers(expr, condition = function(cnd) { . { . .__handler_frame__. <- TRUE . .__setup_frame__. <- frame . if (inherits(cnd, "message")) { . except <- c("warning", "error") . } . else if (inherits(cnd, "warning")) { . except <- "error" . } . else { . except <- "" . } . } . while (!is_null(cnd)) { . if (inherits(cnd, "error")) { . out <- handlers[[1L]](cnd) . if (!inherits(out, "rlang_zap")) . throw(out) . } . inherit <- .subset2(.subset2(cnd, "rlang"), "inherit") . if (is_false(inherit)) { . return() . } . cnd <- .subset2(cnd, "parent") . } . }), stackOverflowError = handlers[[1L]]), tryCatchList(expr, . classes, parentenv, handlers), tryCatchOne(expr, . names, parentenv, handlers[[1L]]), doTryCatch(return(expr), . name, parentenv, handler), withCallingHandlers(expr, . condition = function(cnd) { . { . .__handler_frame__. <- TRUE . .__setup_frame__. <- frame . if (inherits(cnd, "message")) { . except <- c("warning", "error") . } . else if (inherits(cnd, "warning")) { . except <- "error" . } . else { . except <- "" . } . } . while (!is_null(cnd)) { . if (inherits(cnd, "error")) { . out <- handlers[[1L]](cnd) . if (!inherits(out, "rlang_zap")) . throw(out) . } . inherit <- .subset2(.subset2(cnd, "rlang"), . "inherit") . if (is_false(inherit)) { . return() . } . cnd <- .subset2(cnd, "parent") . } . }), f(l = layers[[i]], d = data[[i]]), l$compute_aesthetics(d, . plot), compute_aesthetics(..., self = self), check_aesthetics(evaled, . n), cli::cli_abort(c("Aesthetics must be either length 1 or the same as the data ({n})", . x = "Fix the following mappings: {.col {names(which(!good))}}")), . rlang::abort(message, ..., call = call, use_cli_format = TRUE, . .frame = .frame)), parent = c(0L, 1L, 2L, 3L, 4L, . 5L, 6L, 7L, 6L, 9L, 10L, 4L, 12L, 13L, 14L, 15L, 16L, 17L, . 18L, 13L, 13L, 13L, 22L, 23L, 24L, 25L, 26L, 26L, 25L, 29L, . 30L, 31L, 32L, 33L, 34L, 31L, 31L, 31L, 38L, 38L, 40L, 38L, . 38L, 43L, 43L, 45L, 46L, 47L, 48L, 49L, 50L, 47L, 46L, 53L, . 54L, 55L, 56L, 57L), visible = c(TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, . TRUE, TRUE, FALSE, FALSE), namespace = c("IRkernel", NA, . "IRkernel", NA, "base", "base", "base", "base", "base", "base", . "base", "evaluate", "evaluate", "evaluate", "base", "base", . "base", "base", "base", "base", "base", "IRkernel", "IRkernel", . "IRdisplay", "IRdisplay", "base", "base", "base", "base", . "base", "IRdisplay", "base", "base", "base", "base", "base", . NA, "repr", "base", "utils", "base", "base", "ggplot2", "ggplot2", . "ggplot2", "ggplot2", "rlang", "base", "base", "base", "base", . "base", "ggplot2", NA, "ggplot2", "ggplot2", "cli", "rlang" . ), scope = c("::", NA, "local", NA, "::", "local", "local", . "local", "local", "local", "local", "::", ":::", "local", . "::", "::", "local", "local", "local", "::", "::", "local", . ":::", "::", ":::", "::", "::", "::", "::", "::", "local", . "::", "local", "local", "local", "::", NA, ":::", "::", "::", . "::", "::", ":::", "::", ":::", ":::", "::", "::", "local", . "local", "local", "::", "local", NA, "local", ":::", "::", . "::"), error_frame = c(FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, . FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE)), row.names = c(NA, . -58L), version = 2L, class = c("rlang_trace", "rlib_trace", . "tbl", "data.frame")), parent = NULL, body = c(x = "Fix the following mappings: `x`"), . rlang = list(inherit = TRUE), call = check_aesthetics(evaled, . n), use_cli_format = TRUE), class = c("rlang_error", . "error", "condition"))) 31. handlers[[1L]](cnd) 32. cli::cli_abort(c("Problem while {step}.", i = "Error occurred in the {ordinal(i)} layer."), . call = layers[[i]]$constructor, parent = cnd) 33. rlang::abort(message, ..., call = call, use_cli_format = TRUE, . .frame = .frame) 34. signal_abort(cnd, .file)
Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport" Error in pushViewport(viewport(x = 1, y = 1, width = 0.9, height = 0.9, : could not find function "pushViewport"
# Dimplot1
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(sub.epidermis2, reduction = "umap",label = FALSE,group.by = "seurat_clusters", pt.size = 1.5)
# Dimplot2
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(sub.epidermis2, reduction = "umap",label = FALSE,group.by = "celltype.anno", cols = c("#360167", "#59adf6", "#Ff4061"),pt.size = 1.5)
# Dimplot3
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(sub.epidermis2, reduction = "umap",label = FALSE,group.by = "time.anno", cols = c("red", "grey", "blue"),pt.size = 1.5)
sessionInfo()
R version 4.3.1 (2023-06-16) Platform: x86_64-conda-linux-gnu (64-bit) Running under: CentOS Stream 8 Matrix products: default BLAS/LAPACK: /hpc/group/pbenfeylab/mz187/miniconda3/envs/mon3/lib/libopenblasp-r0.3.23.so; LAPACK version 3.11.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C time zone: America/New_York tzcode source: system (glibc) attached base packages: [1] stats4 stats graphics grDevices utils datasets methods [8] base other attached packages: [1] RColorBrewer_1.1-3 ggplot2_3.4.3 [3] monocle3_1.0.0 SingleCellExperiment_1.22.0 [5] SummarizedExperiment_1.30.2 GenomicRanges_1.52.0 [7] GenomeInfoDb_1.36.1 IRanges_2.34.1 [9] S4Vectors_0.38.1 MatrixGenerics_1.12.2 [11] matrixStats_1.0.0 Biobase_2.60.0 [13] BiocGenerics_0.46.0 SeuratWrappers_0.3.1 [15] SeuratObject_4.1.3 Seurat_4.3.0.1 loaded via a namespace (and not attached): [1] jsonlite_1.8.7 magrittr_2.0.3 spatstat.utils_3.0-3 [4] farver_2.1.1 zlibbioc_1.46.0 vctrs_0.6.3 [7] ROCR_1.0-11 spatstat.explore_3.2-1 RCurl_1.98-1.12 [10] base64enc_0.1-3 S4Arrays_1.0.4 htmltools_0.5.6 [13] sctransform_0.3.5 parallelly_1.36.0 KernSmooth_2.23-22 [16] htmlwidgets_1.6.2 ica_1.0-3 plyr_1.8.8 [19] plotly_4.10.2 zoo_1.8-12 uuid_1.1-1 [22] igraph_1.5.1 mime_0.12 lifecycle_1.0.3 [25] pkgconfig_2.0.3 rsvd_1.0.5 Matrix_1.6-1 [28] R6_2.5.1 fastmap_1.1.1 GenomeInfoDbData_1.2.10 [31] fitdistrplus_1.1-11 future_1.33.0 shiny_1.7.5 [34] digest_0.6.33 colorspace_2.1-0 patchwork_1.1.3 [37] tensor_1.5 irlba_2.3.5.1 labeling_0.4.2 [40] progressr_0.14.0 fansi_1.0.4 spatstat.sparse_3.0-2 [43] httr_1.4.7 polyclip_1.10-4 abind_1.4-5 [46] compiler_4.3.1 proxy_0.4-27 remotes_2.4.2.1 [49] withr_2.5.0 speedglm_0.3-5 DBI_1.1.3 [52] viridis_0.6.4 biglm_0.9-2.1 R.utils_2.12.2 [55] MASS_7.3-60 DelayedArray_0.26.6 tools_4.3.1 [58] lmtest_0.9-40 httpuv_1.6.11 future.apply_1.11.0 [61] goftest_1.2-3 R.oo_1.25.0 glue_1.6.2 [64] nlme_3.1-163 promises_1.2.1 grid_4.3.1 [67] pbdZMQ_0.3-9 Rtsne_0.16 cluster_2.1.4 [70] reshape2_1.4.4 generics_0.1.3 gtable_0.3.4 [73] spatstat.data_3.0-1 R.methodsS3_1.8.2 tidyr_1.3.0 [76] data.table_1.14.8 XVector_0.40.0 sp_2.0-0 [79] utf8_1.2.3 spatstat.geom_3.2-4 RcppAnnoy_0.0.21 [82] ggrepel_0.9.3 RANN_2.6.1 pillar_1.9.0 [85] stringr_1.5.0 IRdisplay_1.1 later_1.3.1 [88] splines_4.3.1 dplyr_1.1.2 lattice_0.21-8 [91] survival_3.5-7 deldir_1.0-9 tidyselect_1.2.0 [94] miniUI_0.1.1.1 pbapply_1.7-2 gridExtra_2.3 [97] scattermore_1.2 stringi_1.7.12 lazyeval_0.2.2 [100] evaluate_0.21 codetools_0.2-19 tibble_3.2.1 [103] BiocManager_1.30.22 cli_3.6.1 uwot_0.1.16 [106] IRkernel_1.3.2.9000 xtable_1.8-4 reticulate_1.31 [109] repr_1.1.6 munsell_0.5.0 Rcpp_1.0.11 [112] globals_0.16.2 spatstat.random_3.1-5 png_0.1-8 [115] parallel_4.3.1 ellipsis_0.3.2 assertthat_0.2.1 [118] bitops_1.0-7 listenv_0.9.0 viridisLite_0.4.2 [121] scales_1.2.1 ggridges_0.5.4 leiden_0.4.3 [124] purrr_1.0.2 crayon_1.5.2 rlang_1.1.1 [127] cowplot_1.1.1